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Using discrete calculus, we derive the missing stress-geometry equation for rigid granular mate- 
rials in two dimensions, in the mean-field approximation. We show that (i) the equation imposes 
granularity in a literal sense; (ii) we recover anisotropic elasticity, without reference to strain, and 
(iii) the packing fabric plays an essential role. 
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INTRODUCTION 

Despite a century of study, the macroscopic behavior 
of quasistatic granular materials remains poorly under- 
stood We still lack a fundamental system of contin- 
uum equations, comparable to the Navier-Stokes equa- 
tions for a Newtonian fluid Q . Experiments and simula- 
tions indicate that stress distribution within a granular 
solid depends on the packing's preparation history Q. 
The latter is known to induce anisotropy in the statis- 
tics of grain arrangement, known as the packing fabric 
[3, Q- Hence the fabric, characterized most simply by 
a 2nd-order symmetric tensor, may be the crucial inter- 
nal variable needed to close macroscopic equations for 
quasistatic granular materials. 

Using tools of discrete calculus, in this work we de- 
rive one of the missing continuum equations for two- 
dimensional granular materials directly from the grain 
scale, in the mean-field approximation. The stress- 
geometry equation thus derived relates the stress tensor 
to the fabric. It unambiguously shows that, subject to 
boundary conditions on stress, frictional granular mate- 
rials are described by a version of anisotropic elasticity, 
as previously suggested 0, [1| . 

The problem is most easily posed for isostatic packings 
of perfectly frictional, rigid grains in the absence of grav- 
ity; later we will relax these assumptions. We consider 
a frictional granular material in the plane, composed of 
Nrq convex grains in mechanical equilibrium, touching 
at Nc point contacts [f| . Mechanical equilibrium requires 
that 
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where f° is the contact force exerted on grain g at contact 
c, r c is the position of c, r 9 is position of the center-of- 
mass of g, and C 9 is the set of contacts belonging to g. 
Here the cross-product is defined &saxb = aib = 
a t Eijbj with £12 = -£21 = 1, Ell = £22 = 0. 

Newton's laws (Q} give 3Nrg scalar constraints on the 
2Nc = Nrgz degrees-of-freedom in the contact forces, 
defining the contact number z. When z = 3, the packing 
is isostatic: given the positions and orientations of the 



grains and the external loading, Newton's laws can be 
solved for the contact forces @, HI . 

The macroscopic object of interest is the stress tensor 
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where G = G(r) is a set of grains centered on the point 
r, occupying the area Ac- 

In principle, the microscopic isostatic contact force so- 
lution can be coarse-grained to produce the macroscopic 
stress tensor. However, this is both computationally and 
analytically intractable. It would be preferable to de- 
termine the macroscopic <x by the solution to contin- 
uum equations. Mechanical equilibrium requires that the 
stress tensor satisfies 



= V • <T, 



er = cr 
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however these 3 equations are insufficient to determine 
the 4 components of dr. In passing from the grain scale 
to the continuum, one macroscopic equation has gone 
missing: the stress-geometry equation [7, 11, 12| . 

To find the hidden equation, we make essential use of 
the voids in between the grains which, in two dimen- 
sions, are uniquely associated with closed loops of grains 
[T3 . [l3| . We will show that the stress-geometry equation 
bears a simple physical interpretation: the voids cannot 
carry any stress. 

A continuum equation can only apply when there is 
a large separation between microscopic and macroscopic 
length scales. The former, denoted by £, is given by the 
grain scale; the latter, denoted, by L, arises from macro- 
scopic boundary conditions: the domain size, or the size 
of the region to which a force is applied. Throughout we 
assume that 8 = £/L <C 1. 



STRESS POTENTIALS 

Here we present a mean-field derivation, omitting de- 
tails which will be presented elsewhere. We define po- 
tentials p and ip such that contact forces f° and torques 
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are written as 
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where I' = £'(c) (£ = £(c)) is the loop to the right of (to 
the left of) the oriented contact c (see Figure [Tb)) , and r l 
is the center of loop I, defined below. Writing the contact 
forces and torques in this way, force and torque balance 
are identically satisfied. Conversely, the latter equations 
are precisely the integrability conditions needed to write 
((!]) and ([5]). The potentials are unique up to an irrele- 
vant additive constant. Equations (Q| and §5§ were first 
written down by Satake [14| . The formulation which uses 
p but not <p was considered by Ball and Blumenfeld [l2| ■ 
The potentials are not independent, since the torques 
computed from p must equal those computed from (p. 
Writing t c t = r c — r e , this imposes Nq constraints 
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Since all other constraints have been satisfied, (j6|) are 
the discrete stress-geometry equations in this formula- 
tion. Our goal is to rewrite these equations in such a 
way that a continuum limit may be taken. Continuum 
expressions are naturally related to sums of discrete ex- 
pressions around closed contours ll|. Since, by Eulcr's 
formula, Nq — Nl + Nrg — 1, with Nl the number of 
loops, it is natural to sum these equations around grains 
and loops to form an equivalent set of constraints that 
are more easily interpreted as continuum equations. 
For example, summing ([6]) around a grain, we find 
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<f drx p= -A 9 (V ■ p) s , (7) 
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where s g circulates anticlockwise around the grain (Fig- 
ure llb[) , and A 9 is the area of the polygon formed by the 
t\ vectors around g. The latter two equalities in ([7]) are 
definitions in discrete calculus [HI]. (V • p) 9 converges 
to its continuum counterpart in the following sense: a 
smooth function p(r) can always be defined such that 
(V • p) 9 - V ■ p(r 9 ) = f 9 : Vp + v^490(VVp), with f 9 
a fluctuating O(l) fabric tensor with zero average [lij ]. 
If A 9 ~ £ 2 and V ~ 1/X, then when averaged over TV 
grains, the relative error in (V • p) 9 goes to zero with 
1/VN and S = ijh. 

By similar manipulations, we obtain a discrete calculus 
expression for the sum of ([6]) around loops. We find 

= (V • p)«, (8) 
= (A^ + (V ■ ((Vp) x P ))' + (V(r x p)f , (9) 

which, together, are the exact discrete calculus reformu- 
lation of ([6]). To establish the relationship between <x, 99, 
and p, we introduce auxiliary variables 
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which, we will see, are the values of the Airy stress func- 
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tion. If we sum (|T0|) around a loop, we see that 
* 1 \- 



cec e 



provided r £ = -p- J2 c ec l rC ' w ith z l the number of con- 
tacts around a loop. Hence (p is nothing but a loop aver- 
age of tj). Again summing (|10[) around a loop, but with 
coefficients l c , we find 



(g'r 1 ■ (v x i,y 



where 
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is another fabric tensor. The second equality in (|13j) in- 
dicates that g converges to Vr = S in the continuum, 
with 6 the identity tensor. Finally, the stress tensor & 
can be written 



(v x P y 



(14) 



Hence, given the values of ijj one can determine tp, p, and 
a. These relations indicate that stress corresponds to 
curvature of the Airy stress function ip. 

The definition (|10[) has a simple geometric interpreta- 
tion (Figure fTaj). In (r,ip) space the plane with normal 
(e ■ p e ,+l) which passes through (r e ,Lp £ ) is described 
by the equation = (r — r , ip — ip e ) ■ (i ■ p e ,+l) = 
(r — r e ) x p l + ip — (p l . Recalling that t\ = r c — r e , 
the definition (fT0|) says that for each loop, we create a 
facet of a surface which passes through (r c , ipf) for each 
contact. The consistency equations ((U are equivalent to 



continuity of ip% at a contact, 



ip^ c y Continuity 



of ip across a contact is equivalent to continuity of the 
surface at that contact. This defines a piecewise linear 
surface with holes at each grain. The latter can always 
be filled in smoothly. Hence, in the continuum, ip is a 
continuous function. 

The introduction of ip already indicates the physical 
meaning of the stress-geometry equation. Indeed, given 
a smooth function ip(r) in the continuum, the necessary 
and sufficient condition that ip c = ip(r c ) yields a discrete 
Airy stress function satisfying (|10p for some p and ip is 
that ip varies linearly across voids. Since stress corre- 
sponds to curvature of ip, this is equivalent to requiring 
that stress is concentrated on the grains. 



CONTINUUM LIMIT 

To obtain continuum equations, we define area- 
weighted averaging operators 

(P)(r) = -^-J2 A9p9 > (Q)(r) = j-Y, Ai Q e > (I 5 ) 

G geG L i<£L 
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(a) Discrete ip surface. Each facet corresponds to 
a loop; colors correspond to the number of 
contacts around a loop. 




(b) Geometry around contact c. 



£ V 



iP{r l 



(c) Cross-section of ip surface, in its original ip (solid) and smoothed ip 
(dashed) versions. The contacts c\ and C2 are part of the loop I. 



FIG. 1: (a) Discrete ip surface, (b) Local geometry, (c) Cross-section of ip surface. 



for fields defined on grains and loops, respectively. The 
sets G and L are neighborhoods around r, which must 
become arbitrarily large in the continuum limit. Each 
average can be considered an expectation value over the 
quenched geometry; terms involving products then in- 
volve correlations. We assume that the geometry is ho- 
mogeneous on a mesoscopic scale L' with £ -c L' <C L, so 
that any correlations involving the geometry can be ne- 
glected. Then (rp) = (r)(p) and (p) = (g)^ 1 ■ (V x ip). 

The preceding equations ©, ©, (32]), (33]), and (flij) 
then give continuum equations g = 5 and 

= Vp (16) 

= A(93-V) (17) 

p = Vxtf (18) 

<r = VxVxV- (19) 

Because p = V x ip, the continuum equation V ■ p = 
is identically satisfied. We conclude that the continuum 
stress geometry equation is A(ip — ip) = 0. It remains to 
establish the continuum relation between tp and ip. The 
discrete relation (|TT]) suggests the naive closure (p = ip; 
however, this would give a continuum equation which is 
identically satisfied. To understand why this relation is 
renormalized, we appeal to the geometric interpretation 
of ip. 

The discrete Airy stress function ip describes a contin- 
uous, but patchwork surface, which is alternately flat and 
curved on voids and grains, respectively. In cross-section, 
this surface appears as the solid curve in Figure [Tc] In 
passing to a continuum, we forget which parts of space 



were voids, and which were grains; the effect of averaging 
is to replace the original patchwork surface with a coarse- 
grained surface ip which is not flat on voids, depicted by 
the dashed curve in Figure [Tc] Each loop shrinks to a 
point, and the loop equation (9]) becomes an equation 
valid at the point r l . The equation A(ip — ip) = thus 
relates <p(r e ) = (ip)(r e ) to ip(r e ) = (ip)(r e ) = ip(r e ). 
Crucially, because contact forces are repulsive, the ip sur- 
face has positive mean curvature; hence ip(r^) = ip{r ) 
will systematically deviate from (p . 

We can estimate ip(r e ') by Taylor expansion. By homo- 
geneity, it is reasonable to force ip(r) to go through all ip c . 
Then, since the coarse-grained surface ip(r) is assumed 
smooth, we can Taylor expand ip(r c ) around ip(r e ), and 
compute (fr = -X J2cec e ^A rC ) exactly, introducing fab- 
ric tensors which characterize the local geometry. Here 
we will fit ip(r) to a quadratic polynomial around a loop; 
higher-order terms are suppressed by powers of S <C 1. 
We have 

{p(r) = i>{r l ) + Ar ■ \7iP(r e ) + \{Ar) 2 : WWip(r e ) (20) 
with Ar = r — r and hence 

cp e = ^(r l ) + \F l : VV^(r £ ), (21) 
defining a fabric tensor 

= Tt E ^ ( 22 ) 
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Denning its average 



(23) 



we expect <p — ip = ^F : Wip in the continuum and 
hence 



= A (P : VV^) , 



(24) 



which is the continuum mean-field stress-geometry equa- 
tion, to leading order in 6. If desired, it can be written 
in terms of a directly using VV"0 = tr(<r)<5 — <x. 



DISCUSSION 

From (fTTj)) we see that the pressure P = itr(<r) = 
^Atp, and hence ^A(V> — f) is the pressure differential 
caused by granularity. The stress-geometry equation has 
the physical interpretation that no such pressure differ- 
ential exists, and thus imposes granularity in a literal 
sense. 

The fabric tensor F contains information about pack- 
ing inhomogeneity and anisotropy. Its trace is approxi- 
mately £ 2 = tr(.F) ~ {{z — 2)</>(l — xo)) _1 , with </>the area 
fraction and xq the fraction of 'rattlers,' grains which are 
trapped in the packing but do not contribute to mechan- 
ical stability. The dominant geometrical dependence is 
through z — 2, which may vary by a factor of two over a 
packing F also describes anisotropyin the contact 
distribution, which develops under shear j4l.[l7|. In terms 
of the more frequently used Fc(r) — jv^c) Segc l°l C ; 
we have approximately F ~ £ 2 £ ■ Fq ■ £ T , implying that 
F and Fq share principal axes. 

In the simplest isotropic and homogeneous case F(r) = 
^£, 2 S, the stress-geometry equation reduces to the bihar- 
monic equation AA?/> = 0, which is the same equation 
satisfied by the Airy stress function tp in isotropic elas- 
ticity (llj . It is noteworthy that we derive this result 
without reference to strain, Hooke's law, or energy. It 
explains the success of isotropic elasticity in the presence 
of an isotropic fabric Si. 

More generally, a homogeneous but anisotropic fabric 
yields the equation described by ip in anisotropic elastic- 
ity [l9j]. Anisotropy induces stretching and rotation of 
stress contours, as observed in experiments For ex- 
ample, the Green's function to a normal forcing at the 
edge of an infinite half-plane is shown in Figure [2j for 

F oc g 9^ ' wn itc triangular region in the plot 

of P has P < 0; here an instability can be produced if 
the forcing is sufficiently strong. 

Subject to boundary conditions on stress and a homo- 
geneous fabric, the stress-geometry equation thus recov- 
ers anisotropic elasticity. To apply boundary conditions 



on displacements would require an analog of Hooke's law 
for rigid grains, so far absent [20| . 



Extensions 

The result (|24[) was derived assuming rigid, perfectly 
frictional grains at isostaticity, in the absence of grav- 
ity. However, none of these assumptions were essential. 
Finite stiffness of grains makes the geometry dependent 
on stress through the constitutive law at contacts. This 
introduces corrections to the fabric on the order of the 
contact deformation, typically less than 10~ 4 of a grain 
diameter in realistic granular materials 

At finite friction, mechanical equilibrium requires that 
each contact force satisfies the Coulomb friction inequal- 
ity < M/Za^ where and are the tangential 
and normal components of the contact force at c, and 
(if is the microscopic friction coefficient. This can be 
written as the pair of inequalities (G c ■ M±) : <r c > 0, 
where G c = n c n c is a fabric tensor constructed from con- 
tact normals n c , and M± = 8 ± — e. Under the same 
mean-field assumptions as earlier, this yields the pair of 
continuum inequalities 



[G ■ M±) : & > 0. 



(25) 



Summing these inequalities implies G : a > 0, a gener- 
alization of P > 0. In the frictionless limit fit 0, (|25|) 
implies (G ■ e) : & = 0, which states that G and <r share 
principal axes. 

The stress-geometry equation needs to be solved sub- 
ject to the Coulomb inequalities (|25|) [21] . If a prospec- 



tive solution violates one of these inequalities, then fail- 
ure must occur within the material; the location of failed 
regions must be tracked by dynamical evolution of the 
material's preparation, beyond the scope of the present 
theory. 

Isostaticity motivates the existence of a missing con- 
tinuum equation through constraint counting. However, 
the uniqueness of the contact forces at isostaticity was 
never used in the derivation. The present work suggests 
that all the relevant information about hyperstaticity is 
contained in the fabric tensor F and boundary condi- 
tions. 

Finally, because (J3|) is linear in stress, body forces are 
easily added at the continuum level. For example, a par- 
ticular solution to a gravitational body force V • &i = g 
is <ti = g ■ rS. Equation (f24)) applies to the homoge- 
neous stress <r, while the Coulomb inequalities apply to 
the total stress & + ct\ . 



CONCLUSION 

To summarize, in this work we have derived the miss- 
ing stress-geometry equation for two-dimensional fric- 



FIG. 2: Stresses resulting from a normal point forcing at the origin of a semi-infinite half plane, with anisotropic 

fabric. From left: a yy , a xx , a xy , P 



tional granular materials, equation (|24p . in the mean- 
field approximation. Stress transmission is governed by 
an equation equivalent to anisotropic elasticity. Subject 
to boundary conditions on stress and knowledge of the 
fabric tensor F, the present theory allows computation 
of the stresses inside a granular solid, provided the ma- 
terial satisfies everywhere the Coulomb stability inequal- 
ities ([25]) . 

The present work raises many questions: (i) When do 
the mean-field assumptions hold? Theoretical arguments 
for frictionless systems indicate that correlation lengths 



diverge as isostaticity is approached |22j, but observa- 
tions in real systems suggest that correlation lengths are 
smaller than 10 grain diameters [23J. (ii) How large a 
system is needed for the present theory to apply? (iii) 
Can a fabric evolution equation be derived with similar 
tools? Preliminary work suggests that anisotropic cor- 
relation functions are necessary, (iv) What happens in 
3D? First steps toward extending the stress potentials 
have been taken, but to extend the present approach it 
appears necessary to give the loops an orientation with 
an anti-symmetric tensor [ill ]. 
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cussions with N.J. Balmforth, R. Blumcnfcld, O. Dau- 
chot, I. Hewitt, J. McElwaine, F. Radjai, and C. Schoof, 
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